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The functional equivalence of the unstructured grid code FUN3D to the the struc- 
tured grid code LAURA (Langley Aerothermodynamic Upwind Relaxation Algorithm) is 
documented for applications of interest to the Entry, Descent, and Landing (EDL) com- 
munity. Examples from an existing suite of regression tests are used to demonstrate the 
functional equivalence, encompassing various thermochemical models and vehicle config- 
urations. Algorithm modifications required for the node-based unstructured grid code 
(FUN3D) to reproduce functionality of the cell-centered structured code (LAURA) are 
also documented. Challenges associated with computation on tetrahedral grids versus 
computation on structured-grid derived hexahedral systems are discussed. 


I. Introduction 

The simulation of aerothermodynamic environments of entry vehicles within NASA over the past two 
decades have been executed with two independently developed codes. The Langley Aerothermodynamic 
Upwind Relaxation Algorithm (LAURA) 1 was first utilized in 1987 for simulations of a proposed Aeroas- 
sist Flight Experiment (AFE) vehicle 2 and later was extensively used in defining the aerothermodynamic 
database for the Mars Pathfinder. The Data Parallel Line Relaxation Code (DPLR), 3 initially developed at 
the University of Minnesota and later adopted and expanded at NASA Ames, has been utilized in parallel 
for defining vehicle entry environments. Both codes have been used for every NASA mission involving at- 
mospheric entry since Mars Pathfinder, including the Space Shuttle and in the analysis of vehicle programs 
including X-33, X-34, X-37 through to the present Multi-Purpose Crew Vehicle (MPCV). Teams at Langley 
and Ames have benefited over the years by the ability to cross check and correct independently developed 
algorithms and physical models used in the two codes - a highly valued capability in situations where ex- 
trapolation from limited ground-based experimental validation to flight conditions is extremely challenging. 

Both LAURA and DPLR are multi-block structured grid codes. Their ability to define aerothermal 
environments around protuberances, scarfed nozzles with jets firing, and other small scale irregularities on an 
otherwise “smooth” vehicle surface generally requires significant investment of time to generate appropriate 
grids. (Overset grid approaches in DPLR alleviate these challenges but still require careful grid generation 
to accommodate resolution of flow topologies ensuing from local interactions.) Unstructured grid algorithms 
offer greater flexibility in both defining complex geometric features and in adapting the grid to the shocks, 
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expansions, and free shear layers in a hypersonic flow. Consequently, NASA has undertaken a program to 
develop and mature the speed and reliability of unstructured grid algorithm capabilities to a point where 
the analyses currently performed by LAURA and DPLR eventually can be replaced. 

The unstructured grid code FUN3D has been modified to include the physical models used in LAURA 
for aerothermodynamic applications. Indeed, both codes share the exact same modules for thermodynamics, 
transport properties, chemical kinetics, thermal relaxation, and radiation. The “generic-gas path” through 
FUN3D has primarily tested a multi-dimensional reconstruction algorithm to address challenges associated 
with obtaining accurate heating on unstructured grids. 4,5 The generic-gas extension of FUN3D is now 
directed on maturing the test suites and user interfaces to enable it to replace LAURA for aerothermodynamic 
simulations, especially focussed on problems supporting Entry, Descent, and Landing (EDL) applications. 
A critical element of this extension is documentation of the functional equivalence of FUN3D to LAURA for 
EDL applications. The objective of this paper is to report on the progress of these functional equivalence 
tests. It will also describe algorithmic modifications required for a node-based unstructured grid code 
(FUN3D) to reproduce functionality of a cell-centered structured code (LAURA). Finally, it will briefly 
discuss new functional capabilities not currently in LAURA that are planned to enhance FUN3D’s utility 
for EDL applications. 


II. FUN3D and LAURA 

FUN3D is a node based, fully unstructured, finite volume solver of the Euler and Navier-Stokes equations. 6 
The baseline method uses Least Squares (LS) gradient information to execute second-order accurate inviscid 
flux reconstruction. ' A Green-Gauss formulation of gradients at the nodes is used for multi-dimensional 
reconstruction. Viscous gradients across elements are also computed from a Green-Gauss formulation. A 
suite of modules includes all of the gas physics models in LAURA and Viscous Flow Algorithm for Complex 
Flow Analysis code (VULCAN) 8 for thermodynamics, transport properties, chemical kinetics, and thermal 
relaxation. The baseline inviscid flux reconstruction algorithms in the generic gas path mimic LAURA in that 
they utilize quasi-one-dimensional (edge-based) reconstruction using Roe’s averaging 9 with Yee’s symmetric 
total variation diminishing algorithm 10 adapted for unstructured grids. 5 

LAURA is a structured-grid, finite-volume solver, specialized for hypersonic re-entry physics. 1,11 Key 
elements of LAURA include Roe’s averaging 9 and Yee’s symmetric total variation diminishing (STVD) 
algorithm 10 for the formulation of second-order, inviscid flux. Yee’s STVD formulation has been found to 
be exceptionally robust and Courant-number-independent using first point-implicit and then line-implicit 
relaxation for hypersonic flow simulations. The STVD algorithm uses a non-linear, min-mod function as 
a flux limiter that maintains second-order accuracy away from extrema but can admit limit cycles in the 
convergence process, particularly in the vicinity of captured shocks. This occurrence usually manifests itself 
as a stalling of convergence at a very low error norm, essentially a benign ringing in the solution at a level 
that has no impact on aerothermodynamic quantities. Viscous flux is computed using central differences. 

Surface heating in FUN3D is captured from the integral conservation of flux through the dual control 
volume bounding the surface in the energy conservation law. That is, gradients of temperature and species 
concentration are directly calculated in the interior of the domain but not across the surface (as in LAURA). 
Rather, the net energy flux through the surface boundary must balance the inviscid and viscous flux across 
all other surfaces of the dual. 

Eigenvalue limiting (Harten entropy fix 12 ) is required by both LAURA and FUN3D to enable stable 
shock captures and prevent the formation the carbuncle phenomenon. 13, 14 Implementation of limiting in 
both codes differs due to infrastructure. 5 In LAURA, it is known that surface heating can be augmented by 
application of eigenvalue limiting across the boundary-layer by as much as 20%. This problem is resolved 
by scaling back the magnitude of the limiter in the direction normal to the boundary layer for the linear 
waves. The unstructured formulation of FUN3D provides no information regarding the orientation of the 
dual control volume surfaces relative to the boundary layer or to the vehicle surface. A viscous switch in 
FUN3D rescales the limiter as a function of local cell Reynolds number. The cell Reynolds number uses a 
characteristic length across a face defined as the ratio of cell volume average across the face divided by face 
area. Impact of this limiting will be discussed after reviewing results of the functional equivalence test cases. 
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III. Functional Equivalence Overview 


LAURA already uses a web-based regression test system that automatically checks for any impact of 
software changes to established functionality. 15 FUN3D uses these established tests (reported in Case A 
- Case O below) to confirm that its own functionality reproduces that of LAURA over a suite of test 
cases that touch most permutations of physical models, total energy, and geometric complexity of interest 
to EDL. These tests span perfect-gas, equilibrium, and thermochemical non-equilibrium gas chemistries 
representative of Earth, Martian, and Titan atmospheres. Axisymmetric and three-dimensional, blunt and 
slender configurations are included. In general, agreement with LAURA on identical grids and identical 
physical models is targeted to +/-4% for surface pressure and heating. In any cases where that target is 
missed additional analysis will be performed to explain the differences (i.e. grid convergence associated 
with node-based versus cell based algorithms or differences in the implementation of limiters required by 
the code infrastructure). Given that the purpose of this paper is to demonstrate functional equivalency 
between LAURA and FUN3D, only relative differences in results between the two codes will be shown 
without providing specific details about vehicle geometry or magnitudes of simulated environments. 

Several capabilities in LAURA work flow have been reproduced in FUN3D to encourage users to make 
the transition. 

1. A “shuffle” utility allows solution restarts from disparate physical models involving different numbers 
of governing equations. This capability, now available in FUN3D, utilizes a template in which new 
variables are initialized as needed, old variables are retained, and variables are eliminated if unused 
in the new model set. Examples include starting an 11-species air model from a 7-species air model 
solution or a two-equation turbulence case from a laminar case. 

2. FUN3D has long had an algorithm for identifying lines of nodes off the body to enable line-implicit 
relaxation. This algorithm has been extended for prismatic grid systems spanning the entire shock 
layer to reproduce the quasi-one-dimensional grid adaptation in LAURA. The outer boundary may be 
aligned with the captured bow shock and the distribution of nodes in the boundary layer may target 
a user specified cell Reynolds number. This infrastructure has also enabled implementation of coupled 
radiation (Case M below) and algebraic turbulence models (Cases F and N below). 

3. Ray tracing from any point on the surface in any direction through the shock layer is available in both 
LAURA and FUN3D to support analysis of uncoupled, three-dimensional radiative transport. The 
capability is required in situations where a tangent-slab approximation is not appropriate for analysis 
of radiative energy transfer to the surface. In the general case (not yet implemented) the algorithm 
must also be able to interpolate information from rays back to solution nodes to enable computation 
of coupled radiation - including the radiation source term in the energy equation. 

All solutions are obtained on identical grids with provenance from the quasi-one-dimensional grid adapta- 
tion algorithm used in the LAURA solutions. These grids follow best practices in LAURA. In all cases, mesh 
spacing An at the wall on the forebody puts the cell Reynolds number in the range 0.1 < p w c w An/ p w < 2 
where 

N ce ii — p w c w An/ p w (1) 

and p w , c w , and p w are the density, speed of sound, and viscosity at the wall, respectively. Requirements for 
a grid-converged solution are problem dependent. Following LAURA best practices and based on previous 
grid convergence studies on closely related problems 16-22 we would expect a maximum change in surface 
heating distribution on the forebody of less than 8% (usually lower for attached flow) . 

Risk in the migration from LAURA to FUN3D has been mitigated through extensive re-use of modules 
currently in LAURA. Modest challenges associated with conversion of cell-based boundary conditions to 
node-based in FUN3D and with reformulation of limiters without dependence on structured-grid infrastruc- 
ture are being addressed. The primary risk to full implementation derives from a continued requirement 
for high quality, semi-structured (prismatic) grids to achieve acceptable accuracy in simulation of heating in 
the presence of strong shocks. Offering options in multi-dimensional reconstruction to allow more extensive 
use of tetrahedral grids and developing options for grid alignment with shocks and/or bow shock fitting 
are mitigating this risk in FUN3D. Ultimately, if requirements for initial grid generation quality are not 
substantially relieved it will be more difficult to convince current users of LAURA and DPLR to switch. 
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IV. Functional Equivalence Tests 

A. Case A: Sphere, Perfect Gas 

Figure la compares the FUN3D solution for pressure p (black line) and laminar, convective heating q (red 
line) to the corresponding LAURA solution (bars). Axisymmetric flow conditions are for a standard, 5 
km/s test case, discussed extensively in previous papers. 4,5 The ordinate of the figure uses a log 10 scale so 
that the LAURA bars, scaled to ±4% of the LAURA solution, appear as a constant height in the figure. 
Pressures and heating predicted by both codes are in excellent agreement. A very slight bump in the heating 
distribution approaching the axis is evident in both solutions is an artifact of eigenvalue limiting to suppress 
a carbuncle. 




(a) Case A, sphere (b) Case B, spherically capped cone 

Figure 1: Comparison of surface pressure and convective heating between simulations using LAURA and 
FUN3D for laminar, hypersonic flow of a perfect gas. 

B. Case B: Spherically-Capped Cone, Perfect Gas 


Flow Field Pressure Difference 




z, m 


(a) Percent difference of pressure in symmetry 
plane. 

Figure 2: Additional functional equivalence metrics 
simulation. 


(b) Percent difference of pressure and heating 
along the surface. 

for Case B, a perfect gas, spherically capped cone 


Similar to Case A, Fig. lb again shows good agreement between FUN3D and LAURA for 1 km/s, laminar 
flow of a perfect gas over a spherically-capped cone. Figure 2 provides more detail of functional equivalence 
between the codes. Figure 2a shows that pressure differences are less than 4% across the symmetry plane 
except for slightly higher differences across the captured shock. Figure 2b shows a heating difference as 
large as 5% in the stagnation region and at the sphere-cone junction. These differences are likely associated 
with node-centered (FUN3D) versus cell-centered (LAURA) formulations of the discretization and heating 
is especially sensitive to such differences. 
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C. Case C: Sphere, Air, Thermochemical Nonequilibrium 


This case reprises the conditions of Case A for a 5- 
species air model in thermochemical non-equilibrium (two- 
temperature model). 11 A radiative equilibrium, super- 
catalytic wall boundary condition (all species catalyzed to 
their respective values of mass fraction in the free stream) 
is applied with emissivity of 0.89. A wall temperature up- 
date factor of 0.05 was sufficient to maintain stability. 

Figure 3 compares the FUN3D solution for pressure p 
(black line) and laminar, convective heating q (red line) 
to the corresponding LAURA solution (bars). Recall, the 
ordinate of the figure uses a log 10 scale so that the LAURA 
bars, scaled to ±4% of the LAURA solution, appear as a 
constant height in the figure. Pressures are in excellent 
agreement. The FUN3D solution skirts just above the 4% 
bar relative to LAURA in the stagnation region - a trend 
that repeats in other non-perfect-gas cases to follow. 

D. Case D: Capsule, Axisymmetric, Laminar 
Flow, Perfect Gas 



Figure 3: Comparison of surface pressure and 
convective heating between simulations using 
LAURA and FUN3D for laminar, hypersonic 
flow of 5-species air in thermal non-equilibrium. 


Entry, descent, and landing vehicles often utilize capsule 

shapes composed of a spherical segment with a circular shoulder transitioning to an after body containing the 
payload. Spherical segment shapes tend to be very blunt and they accentuate potential problems associated 
with carbuncle formation simply because the percentage of the shock layer characterized by small eigenvalues 
is increased. Larger regions of near-stagnation conditions are more challenging to compute. 

Figure 4 compares the FUN3D solution for pressure p 
(black line) and laminar, convective heating q (red line) 
to the corresponding LAURA solution (bars). The flat- 
ter profiles of pressure and heating approaching the axis 
(z = 0) are indicative of the broader stagnation region 
across the spherical segment. The FUN3D surface pres- 
sures are with 4% of the LAURA values at every node 
except on the shoulder, and heating lies along the upper 
boundary of the LAURA 4% bar. Heating and pressure 
around the shoulder (approaching z = 2) is offset by larger 
amounts; however, this difference is easily interpreted as 
interpolation errors between the node-based and cell-based 
solutions across a rapidly varying function. 



E. Case E: Capsule, 3D, Laminar Flow, N 2 as Per- 
fect Gas 

This case simulates a laminar, Mach 8 wind tunnel flow 
over a capsule shape in nitrogen modeled as a perfect gas. 

Like the previous case, it involves flow over a capsule-like 
geometry but at angle-of-attack a of 28 degrees. Figure 
5(a) compares the FUN3D solution for pressure p (black 
line) and laminar, convective heating q (red line) to the corresponding LAURA solution (bars) on a cut 
orthogonal to the symmetry place at the axis to the shoulder. Figure 5(b) presents the same relative 
differences along the symmetry plane on the surface. Figure 5(c) provides a contour of relative differences 
on a color scale that emphasizes locations where FUN3D is more than 4% larger than LAURA (dark red) 
or more than 4% lower than LAURA (dark blue). All three plots show excellent agreement with pressure. 
FUN3D tracks slightly higher than 4% of LAURA for heating over acreage on the upper half of Fig. 5(c) - 
the region farthest from the stagnation point near the lower edge of this figure. 


Figure 4: Comparison of surface pressure and 
convective heating between simulations using 
LAURA and FUN3D for laminar, hypersonic 
(Mach 6), perfect-gas flow over a capsule-like 
geometry. 
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(a) Laminar flow. Planar cut through axis, orthogonal (b) Laminar flow. Symmetry plane cut. 

to symmetry plane. 
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(c) Global view of relative differences in pressure (left) 
and heating (right) for laminar flow. 


(d) Global view of relative differences in pressure (left) 
and heating (right) for turbulent flow. 




(e) Turbulent flow. Planar cut through axis, orthogonal (f) Turbulent flow. Symmetry plane cut. 

to symmetry plane. 


Figure 5: Comparison of surface pressure and convective heating between simulations using LAURA and 
FUN3D for hypersonic (Mach 8), perfect-gas flow over a capsule-like geometry at a = 28. 
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F. Case F: Capsule, 3D, Turbulent Flow, N 2 as Perfect Gas 

This case repeats conditions of the previous case E but models turbulent flow with the algebraic Cebeci- 
Smith 23 ’ 24 turbulence model. Figure 5(e) compares the FUN3D solution for pressure p (black line) and 
laminar, convective heating q (red line) to the corresponding LAURA solution (bars) on a cut orthogonal to 
the symmetry place at the axis to the shoulder. Figure 5(f) presents the same relative differences along the 
symmetry plane on the surface. Figure 5(d) provides a contour of relative differences on a color scale that 
emphasizes locations where FUN3D is more than 4% larger than LAURA (dark red) or more than 4% lower 
than LAURA (dark blue). All three plots again show excellent agreement with pressure. FUN3D heating is 
in better agreement with LAURA for the turbulent simulation though it still tends to run at the high end 
of the 4% band. The red dot at the stagnation point of Fig. 5(d) suggests that differences in the eigenvalue 
limiting algorithm have a stronger effect on the heating than elsewhere on the surface. 




(a) temperature 


(b) pressure 



(c) heating 



(d) Contour plot of relative difference in pressure and 
heating between LAURA and FUN3D. 


Figure 6: Comparison of surface pressure and heating between simulations using LAURA and FUN3D for 
hypersonic (6.35 km/s) flow of 5-species air in chemical nonequilibrium over a capsule geometry at a = 28. 


G. Case G: Capsule, 3D, Laminar Flow, 5-Species 

This case (Fig. 6) simulates hypersonic flow (V)*, = 6.35km/s) over a capsule in 5-species air in chemical 
non-equilibrium at angle-of-attack a of 28 degrees. Species mass fractions at the wall are computed assuming 
a super-catalytic boundary condition. A radiative equilibrium wall temperature boundary is assumed. Line 
cuts along the symmetry plane in Fig. 6(a-c) for temperature, pressure, and heating generally show good 
agreement (within 4%) between FUN3D and LAURA. Details of the comparisons at the shoulders shown in 
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each figure indicate a slightly larger prediction of heating by FUN3D repeating the observation first made 
in case D and also evident in cases E and F. Figure 6(d) provides a contour of relative differences on a color 
scale that emphasizes locations where FUN3D is more than 4% larger than LAURA (dark red) or more 
than 4% lower than LAURA (dark blue) indicating excellent agreement with pressure (left side) and good 
agreement with heating (within 4% over most of the vehicle). A red ellipse at the stagnation point is again 
evident, indicating differences in the eigenvalue limiting algorithm as noted previously in case F. 

H. Case H: Axisymmetric, Equilibrium flow by Free-Energy Minimization 

This case reprises the conditions of Case C for a 5-species air model in chemical equilibrium as computed by 
free-energy minimization. 25 A constant wall temperature (T w = 500K) and super-catalytic wall boundary 
condition are applied. Good agreement is evident between LAURA and FUN3D for pressure and heating in 
Fig. 7(a). Comparison of pressure and temperature in the symmetry plane (Fig. 7(b)) indicate small differ- 
ences associated with the shock capturing. Pressure near the outflow boundary in FUN3D is approximately 
4% lower than in LAURA - an explanation for this difference is not yet evident and was not seen earlier in 
Fig. 2(a). The root mean square of differences across the shock layer is 1.53% for pressure and 0.48% for 
temperature. 




(a) Percent difference of pressure and heating along the (b) Percent difference of pressure and temperature in the 

surface. symmetry plane. 

Figure 7: Comparison of pressure, temperature, and heating between simulations using LAURA and FUN3D 
for laminar, hypersonic flow of 5-species air in chemical equilibrium as computed by free-energy minimization. 

I. Case I: Space Shuttle Orbiter, Finite-Catalytic Wall 

LAURA solutions of hypersonic flow over the Shuttle Orbiter have been reported previously. 26-28 Compar- 
isons of FUN3D with LAURA for this application have also been presented. 5 The comparisons are repeated 
here at Mach 23 in a manner consistent with the metric requirements for these functional equivalency tests. 
This case is the first to test a finite-catalytic wall boundary condition 29 and to test significant vortical 
separation to the lee side of a winged vehicle. 

Comparisons of pressure and heating from LAURA and FUN3D on cuts along the symmetry plane, chord 
wise across the wing, and circumferentially across the wing-body are presented in Fig. 8. The upper leg of 
the curve in these figures show conditions on the wind side of the surface cut, where pressures and heating 
are higher than on the lee side of the vehicle. One exception to this trend is in Fig. 8(b) where heating on 
the leading edge of the vertical tail shows relatively higher magnitude. The wind side surface comparisons 
are generally within the 4% band. The lee side comparisons show significant excursions outside the 4% 
band in both pressure and heating. These lee side excursions in heating are relatively small in absolute 
magnitude. They appear to be in regions where previously separated flow washes over a surface and may 
include integrated effects of algorithm differences on a node distribution that is not grid converged. Some 
differences are caused by a modified outflow boundary to prevent reverse flow off the lee side wing trailing 
edge. 
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(a) Symmetry plane cut - pressure. (b) Symmetry plane cut - heating. 




(c) Chordwise cut across wing - pressure. (d) Chordwise cut across wing - heating. 




(e) Circumferential cut across wing-body - pressure. 


(f) Circumferential cut across wing-body - pressure. 


Figure 8: Comparison of surface pressure and heating between simulations using LAURA and FUN3D for 
hypersonic (M = 23), flow of 5-species air in chemical nonequilibrium over a Shuttle Orbiter geometry at 
a = 40. 
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(a) Planar cut through axis, orthogonal to symmetry 
plane - pressure. 



(c) Symmetry plane cut - pressure. 



(b) Planar cut through axis, orthogonal to symmetry 
plane - heating. 



(d) Symmetry plane cut - heating. 




(e) Symmetry plane cut - wake heating detail. 


(f) Heating difference. Top - windside. Bottom - leeside. 


Figure 9: Comparison of surface pressure and heating between simulations using LAURA and FUN3D for 
hypersonic (M = 10, 3.3 km/s), flow of 5-species air in chemical nonequilibrium over a capsule geometry at 


a = 28. 
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J. Case J: Capsule Including Wake Flow 

This case is the first to feature near wake flow behind a capsule shape including effects of gas chemistry at 
Mach 10, Voa = 3.3km/s. Boundary conditions include a super-catalytic wall and radiative equilibrium wall 
temperatures. Figure 9(a-e) compares pressure and heating between LAURA and FUN3D along planar cuts 
that are orthogonal to or aligned with the symmetry plane. The upper leg of the curves indicate conditions 
on the wind side surface and the lower leg of the curves capture the lee side surface conditions. Both pressure 
and heating are within the 4% band on the wind side (recall the bar height is ± 4% of LAURA and the bar 
height is constant using a logarithmic ordinate). There are significant excursions on the lee side; Fig. 9(e) is 
a rescaled version of Fig. 9(d) to highlight these differences. As with the lee-side surface comparisons of the 
previous Shuttle Orbiter case it appears that the integrated effects of algorithm differences in regions where 
the separated flow may not be grid converged are the likely cause of the excursion. Consider, for example, 
the good agreement evident in the lower left hand corner of Figs. 9(c-d) (upper left hand corner of rescaled 
Fig. 9(e)). This good agreement occurs on the lee side surface where the flow remains attached. It is not 
until the flow separates, and the surface is bathed in flow that has previously separated elsewhere, that the 
differences significantly exceed 4%. At this point it is not known if the failure to achieve grid convergence 
is an issue in only LAURA or only in FUN3D or in both codes. While these percent differences are large 
Fig. 9(f) shows that the magnitude of the differences between LAURA and FUN3D are everywhere less than 
1 W/cm 2 on the fore body (top half of figure) and aft body (bottom half of figure). 

K. Case K: C0 2 

The results of Fig. 10 show good comparisons with pressure and heating for a case designed to bring in 
Martian atmosphere properties into the test suite. This single species C0 2 , thermal equilibrium case does 
not significantly stress the test suite except for the variation of heat capacity of a linear, triatomic molecule. 
More complex chemical models are planned relevant to the Martian atmosphere in the future. 




(a) pressure (b) heating 

Figure 10: Comparison of surface pressure and heating between simulations using LAURA and FUN3D for 
Mach 6.2 (2.8 km/s), flow of single species C0 2 over a sphere. 
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L. Case L: Titan Entry, 18-Species, Thermochemical Nonequilibrium 


In contrast to the simplicity of the chemical model of Case K, the present case tests an 18-species thermo- 
chemical nonequilibrium model for the entry of a 70 degree spherically capped cone into the atmosphere 
of Titan. The atmospheric model assumes 0.9708 N 2 and 0.0292 CH 4 by mass in the free stream. Re- 
maining species include: N, C, H, CH 3 , CH 2 , CH, H 2 , C 2 , NH, CN, Ar, H + , N + , Ar + , Nj, and e _ . 
This is the first test in the suite to include charged 
particles in the model. Boundary conditions in- 
clude a super-catalytic wall and radiative equilib- 
rium wall temperatures. The simulation is for a 
trajectory point in the Titan entry of Mach 22.6 
(5.761 km/s). Axisymmetric flow is assumed. Com- 
parisons of pressure and heating in Fig. 11 indicate 
that FUN3D is within the 4% band of LAURA. 


M. Case M: Coupled Radiation, 11-species 
Air, Thermochemical Nonequilibrium 



A simulation including fully coupled radiative heat- 
ing at the FIRE II 30,31 1643 second trajectory point 
(Uoo = 10.84 km/s and = 0.00078 kg/m 3 ) is 
presented in Fig. 12 . 32 Both simulations used the 
HARA radiation code 16 and an 11-species, two- 
temperature thermo chemical non-equilibrium gas 
model. Agreement between the two codes for con- 
vective and radiative heating is generally good. This 
routine uses the line extraction capability described 
in the previous section to transfer radiative flux in- 
formation between FUN3D and HARA assuming 

the tangent-slab approximation. Matching the radiative heating is simply a confirmation that the pro- 
files of species densities and vibrational-electronic temperatures being passed to HARA from LAURA and 
FUN3D are in good agreement. 


Figure 11: Comparison of surface pressure and heating 
between simulations using LAURA and FUN3D for 
Mach 22.6 (5.761 km/s), flow of 18-species nitrogen- 
methane mixture over a 70 degree spherically-capped 
cone. 


N. Case N: Multi-component Diffusion with Stefan-Maxwell Equations 


The tests in this case engage two new capabili- 
ties. The first is multi-component diffusion through 
a sub-iteration of the Stefan Maxwell equations . 33 
The second is imposition of an equilibrium catalytic 
boundary condition where species mole fractions at 
the wall are determined a s a function of wall pres- 
sure, temperature, and elemental mass fractions dif- 
fusing to the wall . 25 Simulated conditions are for 
Mach 25 flow (8 km/s) of 11-species air in ther- 
mochemical nonequilibrium. Comparison of pres- 
sure and heating between LAURA and FUN3D are 
presented in Fig. 13 for laminar flow and turbu- 
lent flow using the Cebeci-Smith algebraic turbu- 
lence model . 23 Pressures are in excellent agreement. 
The heating in both the laminar and turbulent mod- 
els skirt the upper limit of the 4% band, indicating 
possible effects of the eigenvalue limiting routine dif- 
ferences in the stagnation region. 



Figure 12: Comparison of convective and radia- 

tive heating between simulations using LAURA and 
FUN3D for FIRE II case at 1643 sec. Lines are 
FUN3D. Bars are ± 4% of LAURA. 
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(a) Laminar (b) Turbulent (Cebeci-Smith) 

Figure 13: Comparison of surface pressure and heating between simulations using LAURA and FUN3D for 
Mach 25 (8 km/s) flow of 11-species air over a sphere assuming multi-component diffusion using the Stefan 
Maxwell equations. 


O. Case O: Compression Ramp, Two-Equation Turbulence Model 

Hypersonic, turbulent flow over a 36 deg. ramp at Mach 11.3 (1.769 km/s) is modeled using the Wilcox 2006 
k-ui two-equation turbulence model. 34 This case corresponds to Run 54 from CUBRC 35 and was thoroughly 
studied in a recent paper to quantify simulation uncertainty of shock wave - turbulent boundary layer 
interactions. 36 

Comparisons of pressure and heating are presented in Fig. 14. The pressure distributions are in good 
agreement (within 4%) starting from the strong interaction at the leading edge of the plate at the left through 
the abrupt rise of pressure at the ramp and the post interaction at the right. However, heating from FUN3D 
coming off the strong interaction and post-shock interaction is in poor agreement with the LAURA baseline 
and further analysis of the multi-equation turbulence models are required. 

Recall that the intent of these comparisons is to confirm that, given the same set of physical models, both 
codes give the same result within uncertainty introduced by grid convergence. The accuracy of the simulation 
relative to real physics as determined by experimental measurements with associated measurement error is 
not identifiable in this comparison. Indeed, the CUBRC data and other models show a greater extent of 
separation in front of the ramp than shown in Fig. 14. A compressibility correction is required with this 
model to match the extent of separation. 36 

V. Eigenvalue Limiting 

The sensitivity of stagnation region heating on the magnitude of eigenvalue limiting for case B is demon- 
strated in Fig. 15. The ordinate z defines radial distance from the stagnation point for the axisymmetric 
cone. The abcissa q is the heating magnitude for this perfect-gas case. The combination of line style and 
overlayed symbol shape defines the magnitude of the eigenvalue limiter in the following equations. 

A maximum wave speed is defined 


Amoi = max(|[/|, |U|) + c (2) 

where U , V are velocity components normal and tangent to the cell wall and c is the sound speed. Coefficients 
to limit the eigenvalues are defined 

e a = max [0.01, min (1, A a ,lim(l - s(N ce a)))] (3) 

13 of 21 


American Institute of Aeronautics and Astronautics 




(a) Pressure (b) Heating 

Figure 14: Comparison of surface pressure and heating between simulations using LAURA and FUN3D for 
Mach 11.3 (1.769 km/s) flow of air modeled as a perfect gas over a 36 deg. ramp. 


e c = max [0.00001, min (1, A Cl ;,m(l — s(N ce u)))] (4) 

where N ce u is the cell Reynolds number (see Eq. (1)) and A a ,Um and A c ,Um are user defined limiting coeffi- 
cients. These coefficients are identical in the present work and their values are denoted A in the line style 
legend of Fig. 15. The function s(N ce u ) in Eq. (5) is designed to drive the limiting coefficients e a and e c to 
a minimum value if N ce u < A/, i and to its maximum value if N ce u > \l, 2 - 

s{N ceU ) = 1 - max [0, min (1, ( N ceU - l, 2 - Al,i))] 4 (5) 

The bounds on N ce u are denoted by the ordered pair A l hr the symbol legend of Fig. 15. It is assumed that 
the user has some knowledge of how the grid spacing at the wall was defined and is able to set A l,i and A l ,2 
so that the limiting may be substantially removed near the wall but retained in the upper portions of the 
boundary layer. If limiting is scaled back too much a pre-carbuncle reverse flow in the nodes near the wall is 
evident and a dip in the stagnation region heating is observed. The formulation of the limiting is completed 
as defined in Eqs. (6-9). 


A 0 

A, 


C'a'^max 

&c^max 


( 6 ) 

( 7 ) 
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( 8 ) 

( 9 ) 


For any ordered pair of limiting coefficients on cell Reynolds number A l,i and A l ,2 it is noted that the 
effect of increasing limiting coefficient on the maximum wave speed A is to increase stagnation point heating 
- consistent with LAURA behavior. The effect of the limiter away from the stagnation point (larger z) 
begins to diminish mainly because the \U\ is itself larger and not modified by the limiting. (If the body is 
blunter then the magnitude of \U\ remains small over a greater domain and more affect may be expected.) 
A setting with A = 0.6 (the red lines in Fig. 15) is most consistent with default settings in LAURA. 
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Figure 15: Sensitivity of stagnation region heating to numerical parameters. 


The sensitivity of heating to the limiting based on N ce u is much less than observed in LAURA (possibly 
due to the integral formulation of surface heating discussed previously in Sec. II. The first ordered-pair for 
A l = (0,0.1) (square symbols) does not engage for the given grid with N ce u > 0.1. The second ordered-pair 
for Al = (5, 50) (delta symbols) engages over the lower portion of the boundary layer but produces very little 
affect on heating magnitude. It is not until the third ordered-pair is engaged with A l = (50,500) (gradient 
symbols) that limiting is substantially removed across the domain and an abrupt lowering of heating is 
produced consistent with pre-carbuncle behavior. 

In the tests documented herein eigenvalue limiting settings were chosen in an attempt to minimize the 
effects of limiting without allowing pre-carbuncle heating dips in the stagnation region. On average, it is 
observed that FUN3D solutions for heating tend to run on the high side of the 4% band of the LAURA 
baseline for attached flow. Excursions can be larger for separated flow cases I, J, O. The extent to which 
these larger differences are due to differences in the integrated effect of limiting across the separated shear 
layer or simply due to grid convergence issues is not known. More work is planned to better understand 
these issues and ultimately automate the selection of limiting parameters in this Roe/STVD based algorithm 
and to engage other flux function formulations in the generic gas path of FUN3D. 

VI. Run Times 

A representative time comparison between the two codes is provided in Fig. 16 associated with case H. 
The ordinate in these figures shows the convergence criteria - the percent change in stagnation point heating 
after 1000 relaxation steps. The abcissa documents the number of relaxation steps. Each solution started 
from the same cold start (uniform flow) on the same grid and ended when the change in stagnation point 
heating was less than 0.1%. The machine being used is a desktop workstation (server) with two Xeon 5640 
quad-core processors and 48GB of memory. FUN3D and LAURA were compiled with Intel FORTRAN 
compiler version 12.x for mpi. Six processors were used for both the FUN3D and LAURA solutions. The 
FUN3D time, shown in red, is 472 seconds using a constant Courant number of 5. The LAURA solution, 
run with a default Courant number of 10 , shows point-implicit relaxation times in Fig. 16(a) to convergence 
and switching to line- implicit relaxation after 6000 steps in Fig. 16(b). In the first case, LAURA requires 
449 seconds and in the second case LAURA requires only 117 seconds to execute fewer, more efficient and 
more costly line relaxation steps. 

Inherent differences in algorithm infrastructure make comparisons on equal footing difficult. In this 
axisymmetric case with identical grids, LAURA solves only 64x1x30 cells and FUN3D solves 65x2x31 nodes. 
Both codes freeze Jacobian updates for ten relaxation steps. LAURA updates the transport properties every 
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(a) LAURA point-implicit (b) LAURA point- and line-implicit 

Figure 16: Timing comparisons. 


five steps but FUN3D updates them every relaxation step. It is observed that multi-species cases do not 
converge as efficiently as perfect-gas cases in FUN3D; consequently, it is suspected that the current implicit 
formulation of the multi-species surface boundary condition may impact efficiency of each relaxation step. 


VII. Multi-Dimensional Reconstruction 



Thick lines - FUN3D Hex (STVD) 
Thin lines - FUN3D Tet (STVD) 


Thick lines - FUN3D Hex (STVD) 

Thin lines - FUN3D Tet (MultiDim Recon) 


(a) Edge-based reconstruction on a tetrahedral grid. (b) Multi-dimensional reconstruction on tetrahedral 

grid. 

Figure 17: Pressure distribution in the plane of symmetry for uncoupled radiation simulation for Fire II on 
tetrahedral grid compared to hexahedral grid simulation. 

A simulation related the Fire II test case M, without radiation, highlights some of the algorithmic issues 
related to hypersonic simulation on tetrahedral grids. Comparisons in the symmetry plane for pressure 
and temperature distributions made by FUN3D on tetrahedral and hexahedral grids and using edge-based 
reconstruction and multi-dimensional reconstruction are presented in Figs. 17 and 18. In this case, it is 
noted that a multi-dimensional reconstruction algorithm 5 is required when using tetrahedral grids with a 
symmetric, total-variation diminishing (STVD) formulation. In all cases, the hexahedral grid simulation 
provides simulation quality that is closest to LAURA (not shown here). Nevertheless, it is important to 
provide evidence of simulation quality on tetrahedral grids, as shown here, because grid generation and 
adaptation is most readily implemented on mixed element grid systems including tetrahedra. The importance 
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Thick lines - FUN3D Hex (STVD) 

Thin lines - FUN3D Tet (MultiDim Recon) 



Thick lines - FUN3D Hex (STVD) 

Thin lines - FUN3D Tet (MultiDim Recon) 


(a) Translational-rotational temperature. 


(b) Vibrational-electronic temperature 


Figure 18: Temperature distribution in the plane of symmetry for uncoupled radiation simulation for Fire 
II on hexahedral and tetrahedral grids. 


of multi-dimensional reconstruction on heating has been reported previously. 5 The large scale effect on the 
shock layer evident in Fig. 17(a) is highlighted again here to give example of issues motivating future work 
on shock fitting and reconstruction with basis functions discussed in the next section. 


VIII. Future Work 

As of this writing there are four function capabilities in LAURA that are not yet in FUN3D. These 
include: 

1. RCS jet inflow boundary condition - a boundary feeding a nozzle geometry for which plenum pres- 
sure and temperature are specified and inflow conditions through the boundary are calculated. This 
capability is scheduled for completion prior to October 2013. 

2. Coupled Ablation - a boundary condition in which equilibrium steady state ablation may be assumed 
and ablation rate and species densities at the wall are calculated. This capability is scheduled for 
completion prior to October 2013. More comprehensive coupled ablation with material response codes 
are also being worked in LAURA and plans call for such coupling to be introduced to FUN3D in fiscal 
year 2014. 

3. Micro-aerothermodynamic morphing - First introduced in the X-33 program and later used extensively 
to model Shuttle Orbiter thermal protection system damage, 28 the morphing capability enables a con- 
verged flow simulation over a smooth outer mold line to be captured and provide boundary conditions 
for subsequent aerothermal analysis of flow over a new geometry morphed from the original surface. 
Unstructured grids in FUN3D are more adaptable to morphed surfaces but at present there is no utility 
to capture an existing solution and re-run only the domain around the perturbation. 

4. Grid sequencing with moving shock 

New capabilities (not currently in LAURA but under development for FUN3D) include bow-shock-fitting 
and the use of orthonormal basis functions for multi-dimensional reconstruction. 


A. Bow-shock-fitting 

A shock fitting methodology is being implemented in FUN3D. The initial implementation will be restricted 
to fitting the free stream boundary to the bow shock. This will greatly alleviate the adverse effects associated 
with capturing the bow shock on unstructured grids. The node based unstructured grid implementation will 
involve using a steady state shock fitting method similar to that described by of Bonfiglioli 3 7 which uses the 
Rankine-Hugoniot relations to solve for the jump conditions thereby yielding the post shock state variables 
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GENERIC GAS FUN3D 
Axisymetric Blunt Body 
U=5 km see, Pure N2 
Thermally Perfect Gas K 
One Temperature Model 
Chemically Frozen Mot 
Shock Stand Off -= 0 
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(a) Captured shock. 

Figure 19: A comparison of shock captured and shock htted solutions around a spherical blunt body. 

and the shock speed. The shock speed is used to move the free stream boundary grid nodes in psuedo- 
time and the solution of the linear elasticity equations will ultimately be used to move the remaining grid 
nodes such that they remain conformal to the other surfaces and the interior of the computational domain. 
Initial testing of shock fitting boundary conditions and grid motion strategies have utilized a preexisting 
hex-grid-based line-adaptation algorithm to iteratively move the free stream boundary until it fits the bow 
shock shape. Figure 19 presents a comparison of the static temperature contours from a very preliminary 
initial test of a shock fitted solution with a shock captured solution for 5 km/sec flow of chemically frozen, 
thermally perfect N2 over a spherical blunt body. 

B. Reconstruction with segmented basis functions 

As noted earlier, multi-dimensional reconstruction provides significant improvement to simulations of hyper- 
sonic flow on tetrahedral grids. Solution quality of surface heating is still not as good as that obtained using 
well aligned, hexahedral grids. Sensitivity of the solution to the algorithm detecting a shock normal direction 
is observed in the heating. An approach using segmented, orthonormal basis functions is being developed 
to provide more accurate flux reconstruction within an element that includes a shock with the expectation 
that basis functions comprised of discontinuous segments spanning the element may provide some advantage 
when the solution itself is also discontinuous. 

As these basis functions were derived a self-mapping property under multiplication became evident in 
which the product of any two basis functions gi(x)gj(x ) equals a rescaled basis function [xb — a^) -1 / 2 gk{x) 
where Xb and x a are endpoints of the domain. This property greatly simplifies an algorithm to obtain 
the series solutions of non-linear differential equations. A separate paper on this topic is currently under 
review. Some preliminary results follow to highlight how this method may impact reconstruction algorithms 
in FUN3D. 

The equations defining steady, quasi-one-dimensional flow through a nozzle are 


d(pu 2 A 


dpuA 

dx 

-pA) 


dx 

dpuHA 

dx 


= 0 


= pA 


= 0 


1 dA 
A dx 


(10) 

( 11 ) 

(12) 


where primitive, dimensionless variables are density p, velocity u, pressure p , and internal energy e. The 
total enthalpy H is defined H = e + p/p + u 2 / 2. The equation of state for a perfect gas is p = (7 — l)pe 
where the ratio of specific heats for air as a perfect gas is 7 = 1.4. A superscript * indicates a dimensional 
quantity. Pressure and total enthalpy are non-dimensionalized by their respective values in the nozzle 
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plenum; consequently p p i en = 1 and H p i en = 1. Density is non-dimensionalized by p* = P p i en / H* len so 
that Ppi en = 7/(7 — 1). Velocity is non-dimensionalized by H* len . Distance from the nozzle throat x* is 

non-dimensionalized by the distance to the exit plane x* xit . The nozzle area distribution A(x) used herein 
is defined 

A(x) = A throat + - (1 - COs(7rx)) - 1 < X < 1 (13) 

A nozzle schematic with boundary conditions is shown in Fig. 20. Two inflow boundary conditions at 
x = x a = — 1 match total enthalpy and isentropic flow from the plenum. A single outflow boundary 
condition at x — 27 = 1 matches pressure at the exit plane. 


exit 



Pexit 


Figure 20: Schematic of nozzle geometry and boundary conditions. 


Analytic solutions to nozzle flow can be defined as a function of area ratio. 38 The nozzle has Athroat = 
A(0) = .05 and A ex i t = A{ 1) = 1.05. In this geometry, if 1 > p ex it > 0.9994681 then the flow throughout 
the nozzle is subsonic. If 0.9994681 > p ex it > 0.064722579 then the flow is sonic at the throat and a shock 
stands between the throat and the exit. A simulation with p ex u = 0.1 is presented in Fig. 21. The solution 
to Eqs. (10 - 12) is comprised of basis function representations of 

128 

^ ' Ql,n9n(x) 
n—1 
128 

y ] Q2,n9n(x) 

n—1 
128 

y ' Qs,n9n(x) 

n—1 

across the nozzle. The exact solution for pressure is shown in red and for Mach number is shown in blue. The 
shock is captured as a discontinuity at the correct location in the nozzle. The nozzle domain is not explicitly 
discretized and the solution is derived only as a function of basis function components. The segmented 
nature of the basis function solution is retained in the presentation of the solution (black lines) in Fig. 21. 

Extending this approach to systems of equations in 2 or more directions will be studied for the purpose 
of improving simulations on tetrahedral elements - exploiting the ability to capture discontinuities without 
occurrence of Gibbs phenomena. 


p{x)A{x) = 
p(x)u(x)A(x) = 
p(x)E(x)A(x) = 


IX. Summary 

The functional equivalence of the unstructured grid code FUN3D to the the structured grid code LAURA 
(Langley Aerothermodynamic Upwind Relaxation Algorithm) is documented for applications of interest to 
the Entry, Descent, and Landing (EDL) community. Examples from an existing suite of regression tests are 
used to demonstrate the functional equivalence, encompassing various thermochemical models and vehicle 
configurations. Figures are presented in which a ±4% LAURA band is defined for surface pressures and 
heating. In all cases with attached flow FUN3D pressures go through the center of the band - indicating 
excellent agreement with LAURA. In most cases with attached flow FUN3D convective heating is bounded 
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Figure 21: Series solution to supersonic nozzle flow with a standing shock in the nozzle for A exit j A t h roat = 21, 
128 basis function components, and p ex it = 0.1. The exact solution for pressure (red) and exact solution for 
Mach number (blue) compared to series solution (black). 


by the 4% band or skirts the upper limit of the band in the stagnation region. Radiative heating from the 
FUN3D-HARA coupling is in very good agreement with the LAURA -HARA coupling for a Fire II test 
case. Cases with significant amounts of separated flow showed larger functional differences between FUN3D 
and LAURA for both pressure and heating. While some discussion of the role of eigenvalue limiting in 
these differences is provided - especially for the separated flow case - more work is required to isolate if the 
differences are caused by numerical parameters, a failure to achieve grid convergence in the separated zone, 
or both. 

Algorithm modifications required for the node-based unstructured grid code (FUN3D) to reproduce func- 
tionality of the cell-centered structured code (LAURA) are also documented. These algorithm modifications 
substantially involve boundary condition formulation (with and without availability of ghost cells) and im- 
plementation of eigenvalue limiting when orientation of a cell face relative to the boundary layer cannot 
be assumed based on a structured ( i,j,k ) ordering of cells. Challenges associated with computation on 
tetrahedral grids versus computation on structured-grid derived hexahedral systems are discussed and some 
innovative ideas are introduced to deal with these challenges. 
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